﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

using Nature.Numerics;
using Nature.Numerics.Basic;

namespace Nature.Toolkit.GMix
{
    class ShockWaveSoneFunction : ISoneFunction
    {
        #region Private Fields
        CompositionObject m_composition;
        double m_h0, m_t0, m_ro0, m_p0, m_c0;
        double m_fixedCompressionRatio;
        double m_fixedMachNumber;
        double m_fixedTemperature;
        double m_fixedPressure;
        double m_fixedEnthalpy;
        Action<double[], double[]> m_selectedFormulation;
        #endregion


        internal void Configure(CompositionObject composition, Temperature? temperature, Pressure? pressure)
        {
            m_composition = composition;
            m_t0 = composition.BaseTemperature.Value;
            m_p0 = composition.BasePressure.Value;
            m_ro0 = composition.CalcDensity(m_t0, m_p0);
            m_h0 = composition.CalcEnthalpy(m_t0);
            m_c0 = composition.CalcSoundVelocity(m_t0);

            if (composition.CompressionRatio.HasValue) 
            { 
                m_fixedCompressionRatio = composition.CompressionRatio.Value;
                m_selectedFormulation = this.CompressionRatioFormulation;
            }
            else if (composition.MachNumber.HasValue) 
            {
                m_fixedMachNumber = composition.MachNumber.Value;
                m_selectedFormulation = this.MachNumberFormulation;
            }
            else if (pressure.HasValue && !temperature.HasValue) 
            {
                m_fixedPressure = pressure.Value;
                m_selectedFormulation = this.FixedPressureFormulation;
            }
            else if (temperature.HasValue && !pressure.HasValue)
            {
                m_fixedTemperature = temperature.Value;
                m_fixedEnthalpy = composition.CalcEnthalpy(m_fixedTemperature);
                m_selectedFormulation = this.FixedTemperatureFormulation;
            }
            else
            {
                throw new ArgumentException();
            }
        }

        #region ISoneFunction Members

        public int Nva
        {
            get { return 4; }
        }

        public void Compute(double[] x, double[] y)
        {
            m_selectedFormulation(x, y);
        }

        private void MachNumberFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;

            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_composition.CalcEnthalpy(xv.Temperature);
            double ro = m_composition.CalcDensity(t, p);
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;

            y[0] = (m_h0 - h) + 0.5 * Constants.Rgas * (t / p + m_t0 / m_p0) * (p - m_p0);
            y[1] = j * j - (p - m_p0) / (1.0 / m_ro0 - 1.0 / ro);
            y[2] = cr - ro / m_ro0;
            y[3] = m - m_fixedMachNumber;
        }


        private void CompressionRatioFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;

            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_composition.CalcEnthalpy(xv.Temperature);
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;
            double w = m_composition.MolarMass;

            double fcr = m_fixedCompressionRatio;

            y[0] = (m_h0 - h) / w + 0.5 * m_p0 / m_ro0 * (fcr * t / m_t0 - 1.0) * (fcr + 1.0) / fcr;
            y[1] = m - 1.0 / m_ro0 / m_c0 *
                Math.Sqrt(
                    fcr * m_p0 * m_ro0 * (fcr * t / m_t0 - 1.0) / (fcr - 1.0)
                );
            y[2] = p - m_p0 * fcr * t / m_t0;
            y[3] = cr - fcr;            
        }

        private void FixedTemperatureFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;

            double fxt = m_fixedTemperature;
            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_fixedEnthalpy;
            double ro = m_composition.CalcDensity(fxt, p);
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;
            double w = m_composition.MolarMass;

            y[0] = (m_h0 - h) / w + 0.5 * m_p0 / m_ro0 * (cr * fxt / m_t0 - 1.0) * (cr + 1.0) / cr;
            y[1] = m - 1.0 / m_ro0 / m_c0 *
                Math.Sqrt(
                    cr * m_p0 * m_ro0 * (cr * fxt / m_t0 - 1.0) / (cr - 1.0)
                );
            y[2] = p - m_p0 * cr * fxt / m_t0;
            y[3] = t - fxt;
        }


        private void FixedPressureFormulation(double[] x, double[] y)
        {
            ShockCompressionVector xv = (ShockCompressionVector)x;

            double fxp = m_fixedPressure;

            double t = xv.Temperature;
            double p = xv.Pressure;
            double m = xv.MachNumber;
            double cr = xv.CompressionRatio;
            double h = m_composition.CalcEnthalpy(xv.Temperature);
            double ro = m_composition.CalcDensity(t, fxp);
            double d = xv.MachNumber * m_c0;
            double j = m * m_ro0 * m_c0;
            double w = m_composition.MolarMass;
            

            y[0] = (m_h0 - h) + 0.5 * ( fxp - m_p0) * Constants.Rgas * ( t / fxp + m_t0 / m_p0);
            y[1] = m - 1.0 / m_ro0 / m_c0 *
                Math.Sqrt( 
                        (fxp - m_p0) / (1.0 / m_ro0 - 1.0 / ro)
                    );
            y[2] = cr - ro / m_ro0;
            y[3] = p - fxp;
        }

        #endregion
    }

    class ShockWaveSoneSolver : SoneNewton<ShockWaveSoneFunction>
    {
        internal ShockWaveSoneSolver()
            : base(new ShockWaveSoneFunction())
        {          
            SetAbsoluteConvergenceRadius(
                    ShockCompressionVector.Create(100.0, 10.0 * Constants.Atmosphere, 10.0, 0.5 )
                );
            SetRelativeConvergenceRadius(
                    ShockCompressionVector.Create(1.0e-1, 0.5, 5.0, 0.5)
                );

            SetAbsoluteTolerance(0.0);
            SetRelativeTolerance(1.0e-8);
        }


        
    }
}
